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Abstract 



\ A recent experiment [Deng et al., Nature 398, 218(1999)] demonstrated four- 

c/5 \ wave mixing of matter wavepackets created from a Bose-Einstein condensate. 

The experiment utilized light pulses to create two high-momentum wavepack- 
ets via Bragg diffraction from a stationary Bose-Einstein condensate. The 
high-momentum components and the initial low momentum condensate inter- 
^ , act to form a new momentum component due to the nonlinear self-interaction 

O \ of the bosonic atoms. We develop a three-dimensional quantum mechanical 

description, based on the slowly-varying-envelope approximation, for four- 
wave mixing in Bose-Einstein condensates using the time-dependent Gross- 
^ , Pitaevskii equation. We apply this description to describe the experimen- 

0^ \ tal observations and to make predictions. We examine the role of phase- 

modulation, momentum and energy conservation (i.e., phase-matching), and 
CN ■ particle number conservation in four-wave mixing of matter waves, and de- 

■ velop simple models for understanding our numerical results. 
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I. INTRODUCTION 



Nonlinear optics has been made possible by the nonlinear nature of the interaction be- 
tween light and matter and by the development of intense light sources that can probe the 
nonlinear regime of this interaction. Nonlinear optical processes include three- and four-wave 
mixing (4WM) processes (e.g., second harmonic generation and third harmonic generation). 
In 4WM three waves (or light pulses) mix to produce a fourth. In this paper we detail 
our studies of 4WM of coherent matter waves. Trippenbach et al. proposed a 4WM ex- 
periment using three colliding Bose-Einstein condensate (BEC) wavepackets with different 
momenta. Deng et al. successfully demonstrated 4WM in an experiment with three BEC 
wavepackets, which interact in a nonlinear manner to make a fourth BEC wavepacket. Here 
we greatly elaborate on and further develop the theory and describe numerical simulations 
of the 4WM output that agree well with the experimental measurements of 0. 

The experimental study of nonlinear atom optics is made possible by the advent of 
Bose-Einstein condensation of dilute atomic gases and the atom "laser" [Q, a source 
of coherent matter-waves analogous to the output of optical lasers. A set of optical light 
pulses incident on a parent condensate with momentum Pi = can, by Bragg scattering [Q, 
create two new daughter BEC wavepackets with momenta P2 and P3. Four- wave mixing in 
a single spin-component condensate occurs as a result of the nonlinear self-interaction term 
in the Hamiltonian for a BEC when three such BEC wavepackets with momenta Pi, P2, and 
P3 collide and interact. The nonlinear self-interaction can generate a new BEC wavepacket 
with a new momentum P4 = Pi — P2 -|- P3. 

The possibility of nonlinear effects in atom optics has been long recognized . Goldstein 
et al. P proposed that phase conjugation of matter waves should be possible in analogy to 
this phenomenon in nonlinear optics, including the case of multiple spin-component conden- 
sates 10. They considered the case where a "probe" BEC wavepacket interacts with two 
counter-propagating "pump" wavepackets to generate a fourth that is phase conjugate to the 



probe, where the probe is weak and causes negligible depletion of the pump. Law et al. |T0 
also suggested analogies between interactions in multiple spin-component condensates and 
four-wave mixing. Goldstein and Meystre develop a theory of 4WM in multicomponent 
BECs based on an algebraic angular momentum approach to obtain the modes of the cou- 
pled operator equations. Our treatment for a single spin-component condensate is based on 
the time-dependent Gross-Pitaevskii equation (GPE), which has proved to be highly suc- 
cessful in describing the properties of a variety of actual BEC experiments 0. Thus, our 
treatment is for a zero temperature condensate. It also can describe 4WM with or without 
the presence of a trapping potential. 

The nature of 4WM in BEC collisions of matter waves is unlike 4WM for optical 
wavepacket collisions in dispersive media [p!2|-[T^. The nonlinearity in the case of BEC 



is introduced by collisions rather than by interaction with an external medium, and the 
momentum and energy constraints imposed are different in the two cases. The kinetic en- 
ergy of massive particle waves is quadratic in the wavevector of the particles and given by 
(/ik)^/2m, whereas the energy of a photon is linear in the vacuum wavevector of the photon, 
k, and is given by ^c|k|. Moreover, the momentum of massive particle waves is linear in 
the wavevector of the particles and given by ^k, whereas for light in a dispersive medium, 
it is proportional to the product of the frequency of the light, uj = c|k| and the refractive 
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index, n^u), where the refractive index depends upon frequency (and the propagation direc- 
tion in non-isotropic media). Hence, conservation of energy does not in general guarantee 
conservation of momentum in optical 4WM. Clearly, complications involving the properties 
of an additional medium does not arise in the BEC case. In any case, the creation of new 
BEC wavepackets in 4WM is limited to cases when momentum, energy and particle number 
conservation are simultaneously satisfied. 

In this paper we develop a general three-dimensional (3D) description of four-wave mixing 
in single-spin-component Bose-Einstein condensates using a mean-field approach similar 
to the time-dependent GPE, also known as the nonlinear Schrodinger equation [Q. We 
introduce the slowly- varying-envelope approximation (SVEA), a very powerful tool that not 
only gives insight into the nature of 4WM but also gives a set of four coupled equations 
for the four interacting BEC waves that are more computationally tractable for numerical 
simulations of the time-dependent dynamics. Section |I| explains the experimental situation 
we have in mind and develops the basic theoretical methods. Section |IT1| describes the results 
of our numerical calculations and compares these to the NIST experiment 0. Finally, in 
Sec. we present a summary and conclusion. 



II. THEORY OF MATTER- WAVE FOUR- WAVE MIXING 

In this section we describe the theoretical tools used in our study of 4WM of matter waves. 
Section |1I A| reviews how high momentum components of a BEC can be formed using optical 
Bragg pulses to prepare the initial configuration for the "half collision" event. Section [II B| 
specifies the parameters that describe the strength of the various physical effects that play a 
role in 4WM: diffraction, potential energy, nonlinear self-energy, and collisions between the 
different momentum wavepackets. This Section also describes how to transform between ID, 
2D and 3D calculations involving the GPE. This is important because, without the slowly- 
varying-envelope approximation (SVEA) that we introduce below, full 3D calculations are 
too computationally expensive to carry out for the actual experimental conditions. Hence, 



the SVEA must be explicitly checked in 2D against the full CP solution. Section C 



describes the details of the SVEA approximation for 4WM. Then Section |IID| introduces a 



simple estimate for the 4WM output. Finally, Section |IIE| shows how the effect of elastic 
scattering between atoms in different momentum wavepackets can be accounted for. This 
process causes loss of atoms from the wavepackets and lowers the 4WM output. 

Let us consider three BEC wavepackets moving with central momenta Pi, P2, and 
P3. Such moving wavepackets can be created, for example, by optically-induced Bragg 
diffraction of a condensate [Q. If these three wavepackets overlap spatially, the self-energy 
of the atoms can produce matter-wave 4WM, just as the third-order Kerr type nonlinearity 
can produce optical 4WM in nonlinear media. One can imagine a number of scenarios in 
which 4WM can occur in matter-wave interactions. One can consider a "whole collision" 
in which three initially separated BEC wavepackets collide together at the same time, or 
a "half collision" in which the wavepackets are initially formed in the same condensate at 
(nearly) the same time. Although we considered the "whole collision" case in Ref. the 
"half collision" case is easier to realize experimentally using the above-mentioned Bragg 
diffraction technique 0. In what follows, we consider only this configuration, in which the 
three wavepackets initially overlap because they have been created as copies of the initial 
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condensate. These wavepackets have different non- vanishing central momenta and therefore 
they fly apart from one another after they have been created. 

Fig. |I]a shows the basic configuration in momentum space of the wavepackets which we 
consider here. Two daughter condensate wavepackets with momenta P2 and P3 are created 
from a parent condensate with mean momentum Pi = 0. Fig. ^ shows these three momenta 
in the lab frame in which the experiment is carried out at two different times: during the 
early stage of the "half collision" when they still overlap spatially, and at a later time when 
they have spatially separated into four distinct wavepackets. We let P3 lie along the x-axis 
of the coordinate system, and P2 make some angle 6 with respect to the a;-axis. Nonlinear 
4WM creates a fourth wavepacket with momentum P4 = Pi — P2 + P3. We demonstrate 
below in Sec. |1I C| that four-wave mixing of matter waves is only possible if there exists 
a coordinate frame in which the mixing is degenerate, that is, all four P^ values in this 
frame have the same magnitude. Fig. shows the degenerate frame corresponding to a 
moving frame with velocity ^deg = (Pi + P3) / (2m), where m is the atomic mass. The total 
momentum is zero in the degenerate frame, and the wavepackets move in oppositely moving 
pairs. The angle 9' between the vectors P2 and P3 is arbitrary. In the laboratory frame, the 
angle 9 is given by 6' = 9' /2, and the length of the vector P2 is given by IP2I = IP3I cos(6'). 
Fig. |I]b shows a set of different possible values of P2. 



A. Bragg Pulse Creation of High Momentum Components 

We assume that the condensate has only a single spin-component, and that its dynamics 
can be described by the GPE, which is known to provide an excellent account of condensate 
properties 0]: 

'^^^W^ = (^r + V(r, t) + iVt/o|^|2)^(r, t), (1) 

where = -^V^ is the kinetic energy operator, V{r,t) is the external potential imposed 
on the atoms, NUq = A^ ^™°^ is the atom-atom interaction strength that is proportional 
to the s-wave scattering length uq (assumed to be positive), m is the atomic mass, and A'" 
is the total number of atoms. The numerical methods for solving the GPE are described 
below in Sec. |T|. 

First, we use the GPE to obtain the ground state condensate in the trapping potential 
at time t = 0, \l/(r, t = 0). This condensate wavefunction is centered around r = 0, and 
normalized to unity. We assume, as is the case in the NIST experiments 0, that the trapping 
potential V{t, t) is turned off at t = and that the condensate is allowed to evolve under 
the influence of only the mean- field interaction until time ti. This includes the special case 
ti = 0. We could equally well treat the case of leaving the trap on, and we would obtain 
similar results. Eq. (|l|) determines the evolved condensate wavefunction, \l/(r, ti). After 
this period of free evolution, the Bragg pulses are applied to create the wavepackets with 
momenta Pi, P2 and P3. The momentum differences |Pj — Pj\ are much larger than the 
momentum spread of the initial parent BEG wavepacket. The experimental time scale 6t 
for creating these wavepackets is short (~ 70 /is) compared to the time scale on which the 
wavepackets evolve. The state at ^2 = ^1 + provides the initial condition for subsequent 
evolution of these three wavepackets as they undergo nonlinear evolution. 
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The initial state at t2 immediately after the Bragg pulse sequences can be approximated 
in a number of ways. In principle one could set up a set of coupled GPEs for the ground 
and excited atomic state components and explicitly include the effect of coupling the light 
field to the excited electronic state. A simpler approach would be to carry out an adiabatic 
elimination of the excited state and develop an effective light-shift potential in which the 
ground state atoms move. If such approaches are carried out in this case, they show that 
the light acts as a "sudden" perturbation such that each of the wavepackets with central 
momenta Pi, P2 and P3 is to a very good approximation simply a "copy" of the parent 
condensate at t = ti |]T3|. Thus, the initial condition immediately after the application of 
the Bragg pulses can be approximated as being comprised of three BEG wavepackets, 

*(r, t2) = v|/(r, h) ^ exp(^P, ■ r/h), (2) 

i=l 

where fi = Ni/N is the fraction of atoms in wavepacket i, and J2i=i /i = 1 so the norm of 
\l/ remains unity. 

After the formation of the wavepackets with momenta Pi, P2 and P3, the initial wave- 
function in Eq. (|^) evolves, and the wavepackets with the different momenta separate. Dur- 
ing this separation,the nonlinear term in the GPE generates a wavepacket with central 
momentum P4 = Pi — P2 + P3, as long as the constraints discussed in relation to Figs. |l| 
and are satisfied. Energy and momentum are conserved during the wavepacket evolution. 
This can be readily checked by verifying that dE{t)/dt = and dP{t)/dt = 0, where 

E{t) = {^mTr + \uom')Mt)), (3) 
is the energy per particle and 

P{t) = -th{^{t)\V\^it)) , (4) 

is the momentum per particle. We have verified numerically that energy and momentum 
are indeed conserved in our calculations described in Section pl| . 



B. Characteristic Time Scales, and Dimensionless Parameters 

In this subsection we discuss characteristic time scales that can be used to estimate 
the importance of the various effects occurring during the dynamics for a particular set of 
experimental parameters. It is convenient to use the Thomas-Fermi (TF) approximation |^ 
to give quantitative estimates of the size of the condensate and the time scales characterizing 
the dynamics. In the TF approximation, one neglects the kinetic energy operator in the 
time-independent nonlinear Schrodinger equation, 

/i^ = (Tr + \/(r,t) + iVf/o|^n^, (5) 

where n is the chemical potential, to obtain the following analytical expression for the 
wavefunction: |\E'(r)p = ^^^^^ for r such that V^r) < yU and '^{r) = otherwise. The TF 
approximation is valid for sufficiently large numbers of atoms N. It is convenient to define the 
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geometric average of the oscillator frequencies for an asymmetric harmonic potential as = 
{uxUJyiUzy^^. The size of the condensate is then given by the TF radius rxF = •\/2/i/(ma;), 



where the TF approximation to the chemical potential /i is determined by the normalization 
of the wavefunction to unity and is given by /i = | (^^4^) {muj'^)^^^. Hence, the TF 
radius ttf scales with as N^^^. The size of the TF wavepacket in the i = x, y, and z 
directions is TTpii) = {l^/^i)^TF- 

In order to estimate the importance of the various terms in the GPE, we set V = for 
free wavepacket evolution and rewrite Eq. in terms of characteristic time scales tj:,p for 
diffraction, and tjvi for the nonlinear interaction, in the following manner 



'dt 



r^p , , , d\ 1 1^1 



toF dx"^ dy"^ dz^ t^i |^ 



m 



|2 



VI/. (6) 



The diffraction time and the nonlinear interaction time are given by Ibf = ^mr^p/h, t^L = 
(A^f/o|^mP/^)~^; respectively. Here |\l/mP is the maximum value of |\l/(r)p, i.e., |\l/mP = 
|\I^(0)p; hence in the TF approximation, = n/h. The smaller the characteristic time, 
the larger is the corresponding term in the GPE. We also define the collision duration time 
tcoi = {'^i^tf)/v, where v = (P3 — Pi)/m is the initial relative velocity of wavepackets 
1 and 3. Thus, tcoi is the time it takes the wavepackets 1 and 3 to move so that they 
just touch at their TF radii, and therefore no longer overlap. The ratio tcoi/tNL gives an 
indication of the strength of the nonlinearity during the collision. The larger the ratio of 
tcoi/tNL, the stronger the effects of the nonlinearity during the overlap of the wavepackets. 
These characteristic times stand in the ratios tnp '■ tcoi '■ tNL = ^ '■ 2-KrTF ' ial^^ where 
A is the De Broglie wavelength associated with the wavepacket velocity v. Experimental 
condensates with tcoi/tNL ^ 1 can be readily achieved. Thus, the nonlinear term will have 
time to act while the BEG wavepackets remain physically overlapped during a collision. 
Another relevant time scale in the dynamics is the characteristic condensate expansion time, 
texp = In the typical experiments modeled below, t^F ^ ^eip > tcoi > tNL- 

In addition to time scales, there are several natural length scales that are important: the 
size ttf of the condensate, the scale (Afc)~^ of phase variation across the parent condensate 
as it expands and develops a momentum spread h^k due to the mean field potential, and the 
scale {k')~^ of phase variation due to the fast imparted momentum P' = hk', where P' is the 
common magnitude of the momentum for the packets in the degenerate frame (Fig. |l]). These 
stand in the relation {k')~^ -C {Ak)~^ -C ttf- The grid spacings in numerical calculations 
are determined by the necessity to resolve the wavefunction on its fastest scale of variation. 
Thus, using the form of Eq. (0) for \l/ requires a grid smaller than {k')~^. This requirement 
limits practical calculations to 2 dimensions (2D). We will introduce an approximation in 
the next section that allows three dimensional (3D) calculations by eliminating the rapidly 
varying phase factors from the equations to be solved. 

We find it convenient to use reduced dimensionless variables to calculate the dynamics. 
The most commonly used set of reduced dimensionless variables in BEG problems involves 
using "trap units" @]. Here however, except for determining the initial conditions at t = 0, 
the trap potential is turned off, and trap units are not particularly relevant. Since we do 
both 2D and 3D calculations, some care is needed in developing a set of units. The primary 
requirement to simulate 3D experiments with a 2D model is that the relations between the 
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characteristic timescales, top, tcoi and t^L, are as determined by experiment. We have done 
this by scahng the solution of the d-dimensional time-dependent GPE by a (i-dimensional 
volume so that the coefficient of the nonlinear term depends only on the dimension and 



the chemical potential /i. By scaling the condensate wavefunction as \E' = "^/yrj^p, the d- 
dimensional time- dependent GPE for a harmonic potential with frequencies Uj, j = 1 . . . d 
can be written as 

Here ^ is dimensionless for any d, and the known fipp for the 3D problem can be transferred 
to an equivalent time-dependent GPE for a 2D calculation. Furthermore, if we define the 
reduced unit of length, xr, to be xr = ttf, define the unit of time, tji, such that = 
mx\/ {2h), and use the normalization condition: / \^!\'^ d'^r / x'r = 1, we preserve the ratios 
between the most important time scales of the problem. The nonlinear time scale, tNL 
depends only on ^tf and is independent of dimension. The specific relations between the 
3D nonlinear coupling parameter U^^ multiplying (r) \E' (r) in Eq. (0) and U^^ and Uq^ , 

the respective self-energy parameters in ID and 2D are: Uq^ = ^U^^, and U^^ = ^Uq^ . 
These values for [/q insure that the chemical potential ^tf (and all the time scales) are the 
same as in 3D. 



C. Slowly Varying Envelope Approximation 

Let us consider the case when the total wavefunction consists of four wavepackets moving 
with different central momenta Pj = hki, z = 1, . . . , 4. We write the wavefunction as 

4 

^(r, t) = E Mr, t) exp [«(kir - ujit)] , (8) 

i=l 

in order to separate out explicitly the fast oscillating phase factors representing central 
momentum ^kj and kinetic energy Ei = fioji = h^kf/2m. The slowly varying envelopes 
$j(r, t) vary in time and space on much longer scales than the phases. The number of atoms 
in each wavepacket is iVj = jy |$j(r, t)p(i^r, and Ylt=i-^i = N is a. constant. Although 
the slowly varying envelope (^^{v^t = 0) is unpopulated initially, it evolves and becomes 
populated as a result of the 4WM process. If we substitute the expanded form of the 
wavefunction in Eq. (§) into the GPE, collect terms multiplying the same phase factors, 
multiply by the complex conjugate of the appropriate phase factors, and neglect all terms 
that are not phase matched (phase matched terms have stationary phases, do not oscillate, 
and satisfy Eqs. ( p!0| - |TTD below), we obtain a set of coupled equations for the slowly varying 
envelopes $i(r, t): 

6{uJi + uji* — ujj — ujj*) X 
%.(r,t)<|.*.(r,t)<|.,(r,t) , (9) 



^1 + {hkjm) . V + 'j-i-^^' + V{r, t))) $,(r, t) 
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where the delta-functions represent Kronecker delta-functions that are unity when the ar- 
gument vanishes. Mixing between different momentum components can result from the 
nonvanishing nonlinear terms in Eq. which satisfy the phase matching constraints re- 
quired by momentum and energy conservation: 



ki + k,. - kj - kj, = 0, (10) 

^2 _ L.2 



kl + kl -k^-e = 0. (11) 



Each of the indices may take any value between 1 and 4. Eqs. ([T0|) and (|TT]) 

are automatically satisfied in two cases: (a) i = i* = j = j* (all indices are equal), or 
(b) j = i j* = i* (two pairs of equal indices). The corresponding terms describe what 
is called in nonlinear optics cross and self modulation terms respectively. The cross and 
self phase modulation terms do not involve particle exchange between different momentum 
components. In the absence of the trapping potential they modify both amplitude and phase 
of the wavepacket through the mean field interaction. Particle exchange between different 
momentum wavepackets occurs only when all four indices in Eq. (^ are different, and 
conservation of momentum and energy of the atoms participating in the exchange process 
occurs. A set of coupled equations involving wave mixing between the various momentum 
components is therefore obtained. 



The momentum conservation of Eq. (10) implies kj + kj. = kj + kj. = k. It is always 
possible to construct a special reference frame, which we call the degenerate frame, where 
K= 0. Consequently, in this frame kj = — kj. and kj = — kj.. In addition energy conservation 
in Eq. (|Tl|) imposes the condition |kj| = |kj| in the degenerate frame. In this frame all four 
momenta are equal in magnitude and can be divided into two pairs of opposite vectors. 
This explains the use of the conjugated pairs of symbols and (j, j*) in our notation. 

The total number of particles, in all wavepackets, is a conserved quantity. The geometrical 
configuration of the wavepacket momenta in the degenerate frame are illustrated in Fig. ^b. 
In the figure we see two pairs of conjugate wavepackets (1,3) and (2,4). All four momenta are 
equal in magnitude and momenta P'^ and Pg are opposite as are the momenta P2 and P4. 
The angle 6 depicted in the figure is completely arbitrary. However, ^ is not allowed, 
since the wavepackets would no longer be distinguishable. Fig. ^ shows a range of possible 
P2 values for wavepackets in the lab frame that satisfy the phase-matching conditions in 
Eqs. ( [T(]| ) and (|TTp. These conditions only allow IP2I = IP3I cos(6'). 

4WM can be viewed as a process in which one particle is annihilated in each wavepacket 
belonging to an initially populated pair of wavepackets and simultaneously one particle 
is created in each of two wavepackets of another pair, one of which is initially populated 
and the other (wavepacket 4) is initially unpopulated. Hence, using Fig. ^ in the moving 
degenerate frame, 4WM removes one atom from each of the "pump" wavepackets 1 and 
3, and places one atom in the "probe" wavepackets 2 and one atom in the 4WM output 
wavepacket 4. This picture is a consequence of the nature of the nonlinear terms in the four 
SVEA equations. It is this bosonic stimulation of scattering that mimics the stimulated 
emission of photons from an optical nonlinear medium. 

The full SVEA equations for 4WM are explicitly given by: 

[I + {hk,/m) . V + ^(^ + V{r, t))] $i(r, t) = 
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-^A^t/o(|$i|' + 2|<l>2|' + 2|$3|2 + 2|$4n$i - ^NUo<l>4^2<^l , (12) 
(I + {nk,/m) . V + ^(^V^ + l^(r, t))^ $2(r, t) = 

-^Ar[/o(|<l>2|' + 2|<l>i|2 + 2|$3|' + 2|<l>4n$2 - ^iVf/o$:$i$3 , (13) 

I + (/.ka/m) ■ V + + V(r, t)) j $3(r, t) = 

-^iV?7o(|$3r + 2|$i|' + 2|$2|' + 2|<l'4n$3 - ^Arf/o<l>4$t$2 , (14) 

I + (/.k4/m) ■ V + ^(^V^ + r (r, t))^ $4(r, t) = 

-^Ar[/o(|<l>4|' + 2|<l>i|2 + 2|<l>2p + 2|<l>3n<l>4 - ^Art/o$i$;<l>3 . (15) 
n h 

The left hand side of these equations describes the motion of the wavepackets due to their 
kinetic and potential energies. The right hand side describes the effect of the phase matched 
nonlinear interaction terms. The last term on the right hand side of each of the SVEA 
equations is a source term which either creates or destroys atoms in the wavepacket being 
propagated. The other terms on the right hand side of the equations account for the self- 
and cross-phase modulation. These phase modulation terms provide an effective potential 
for each wavepacket that accelerates the atoms in it and modifies its internal momentum 
distribution. 

Before we propagate the SVEA equations, the initial wavefunction of the parent conden- 
sate is determined using the time-dependent GPE. First, the propagation is in imaginary 
time to obtain the initial eigenstate in the presence of the magnetic potential. Then, after 
turning off the magnetic potential, the free evolution in the absence of a trapping potential 
is calculated to provide the initial condition in Eq. (^. This free evolution causes a spatially 
varying phase to develop across the condensate as it expands in the absence of the trapping 
potential. Given the initial condition, the SVEA equations can be used to propagate the 
envelope function of each wavepacket, using the same numerical method used to propagate 
the ordinary time-dependent GPE. 



D. Simple Approximations and Scaling with N 

An estimate of the number of atoms that will be transferred to the 4WM wavepacket 
can be developed as follows. To get the small signal growth at early times, multiply both 
sides of the dynamical equation for the rate of change of $4 , where for simplicity we keep 
only the 4WM term on the right hand side of the equation, 

d^A i 
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by a small time increment 5t to get the growth (5$4 in $4 during 5t: 



^ -i{hf2hf/'^m'^5t ^ -lihhhY"^^ . (17) 

11 iNL 

Here /j = Ni/N is the initial fraction of atoms in wavepacket i, and we assume that $j = 

1 /2 

/j \1/ at early times, because the three wavepackets initially satisfy this relation. Since 
most of the growth takes place in the center of the packets where \1/ is the largest, the factor 
NUq\^\'^ /h is approximated by l/t^L = NUq\'^{0)\'^ /h. Upon squaring this equation, and 
integrating over all space, the total growth in the 4WM output 5/4 is 

Thus, the 4WM signal should grow quadratically at early times. If we take 6t to be the 
total interaction time tcoi defined in Section |11B| , then an estimate of the total 4wm output 
fraction is 

/4 = ^^-/l/2/3(^)'. (19) 

This should be an upper bound on the 4WM output, since the mutual interaction of the 
packets due to the self- and cross-phase modulation terms (the self- and cross-interaction 
energy terms), and their separation from one another when t ~ tcoi, will lower the output. 
Using the TF approximation, = n/h ~ A^^/s ^-^^ ^^^^ _ 2rT^/t> ~ N^^^. Thus, the 

output fraction jf- ~ (^N^/^N'^/^y scales as N^^^. This scaling, which was discussed in 
reference 01, will be checked in our numerical calculations below. 



E. Elastic scattering loss 

Atoms from two different momentum wavepackets can undergo s-wave elastic scattering 



that removes the atoms from the packets and scatters them into An steradians [20|. This 
becomes important when the mean-free-path imfp becomes comparable to or smaller than 
the condensate size, ttf- The mean-free-path is imfp = (c"n)~^, where a = Sna^ is the elastic 
scattering cross section and n is the mean density. Profuse elastic scattering of this type has 
been recently observed [^. This mechanism can also affect the 4WM process since loss of 



atoms from the moving packets reduce the nonlinear source terms in the SVEA equations. 
Although the cloud of elastically scattered atoms can not be simply described by the mean- 
field picture, the loss of atoms from the wavepackets due to this elastic scattering mechanism 
can be described in terms of the SVEA. This is because each momentum component is treated 
separately, and the loss terms due to elastic scattering can be added to the SVEA equations. 

The elastic scattering loss is incorporated by adding loss terms to the right hand side 
of the envelope equations in the form of imaginary potentials that are proportional to the 
density of the "other" momentum component involved in the elastic scattering. The full 



SVEA equations for 4WM, including the effects of elastic scattering loss [20|, are given by: 
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|$l|^$4-^^ — |$2| $4-^^ — |$3r$4 • (23) 

z, 2 2 

There are three elastic scattering loss terms for each SVE momentum component $j arising 
from the interaction of each momentum component with the other three momentum com- 
ponents. The factor of | in the loss terms is due to the fact that these are equations for the 
amplitudes, not the densities. 

The density dependence of the elastic scattering loss terms is identical to that of the 
mean-field interaction terms since both terms are due to elastic scattering. It is of interest to 
compare the strength (size of the coefficient) of the loss term due to elastic scattering with the 
nonlinear term in the GPE. The nonlinear term has a coefficient Uq/H = 4nhao/m, whereas 
the loss term for interaction of packets i and j has a coefficient |fcr = 4nh\ki — kj\aQ/m, 
where v is the relative velocity. The ratio TZ = {^va)/{Uo/h) of loss to mean-field terms for 
packets 1 and 3 in Fig. |l| is 

7^ = 2|kl|ao. (24) 
This ratio is about 0.06 for the NIST 4WM experiment M. 
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III. NUMERICAL SIMULATIONS 



A. Experimental Configuration 

In the NIST experiment 0, the initial sodium F, Mp = 1, —1 condensate is comprised of 
magnetically confined atoms in a TOP (time-orbiting-potential) trap without a discernible 
non-condensed fraction. The trap is adiabatically expanded to reduce the trap frequencies 
in the x, y and z directions to 84, 59 and 42 Hz (the frequency ratios are Ux ■ ujy : ujz = 
1 : 1/v^ : 1/2). After adiabatic expansion, the trap is switched off by removing the 
confining magnetic fields. The condensate freely expands during a delay time ti = 600 fis, 
after which a sequence of two Bragg pulses of 589 nm wavelength creates the two moving 
wavepackets 2 and 3. Each 30 fis Bragg pulse is composed of two linearly polarized laser 
beams detuned from the 3S'i/2, F = 1, Mp = — 1 — 3P3/2, F = 2, Mp = 2 transition by 
about A/2'7r = —2 GHz to suppress spontaneous emission and scattering of the optical 
waves by the atoms. The frequency difference between the two laser beams of a single 
Bragg pulse is chosen to fulfill a first-order Bragg diffraction condition that changes the 
momentum state of the atoms without changing their internal state. The first Bragg pulse 
is composed of two mutually perpendicular laser beams of frequencies Ua and 1^/3 = 1^0 — 50 
kHz, and wavevectors and = fcx and k/s = ky. This pulse sequence causes a fraction /2 
of the BEG atoms to acquire momentum P2 = ^(k^ — k^) = hk{'k + y). A second set of 
Bragg pulses is applied 20 ms after the end of the first Bragg pulse sequence. This pulse is 
composed of two counter-propagating laser beams with frequencies z/q, and I'p = I'a — 100 
kHz, and wavevectors and k^ = A;x and k/3 = — A;x. This pulse sequence causes a fraction 
/s of the BEG atoms to acquire momentum P3 = h{\<.a — k/3) = 2^/cx. Thus, there are 
three initial condensate wavepackets with momenta Pi = 0, P2 and P3 as shown in Fig. |l|. 
The respective wavepacket populations, /i = 1 — /2 — /s, /2, and /s, have a typical ratio 
/i : /2 : /3 = 7 : 3 : 7. 

The number of atoms could be varied between around 3 x 10^ and 3 x 10^. As a typical 



example, we take = 1.5 x 10 atoms in the trap. Taking qq = 2.8 nm ||22|, the nonlinear 
time is t^p = 96.2 /is. The Thomas Fermi radius is rpF = 20.3 fim. Since the separation 
velocity defined in Section [II B| is f = 0.0691 m/s for light of wavelength 589 nm, the physical 



separation time tcoi = ^^-^ = 687 /is in the NIST experiment, and indeed is longer than 
the nonlinear time. The characteristic condensate expansion time, texp = = 1-89 ms for 
a trap with u = 27r-^ s~^. The characteristic diffraction time tp)p = 2mr^p/h = 300 ms 
provides by far the longest time scale in the dynamics. Thus, there is negligible diffraction 
on the time scale of the experiment. 



B. Simulations of the NIST Experiments 

Our solution to the time-dependent GPE uses a standard split-operator fast Fourier 
transform method to propagate an initial state forward in time p3|- The initial state \&(r, t = 
0) of the condensate in the trap is found by iteratively propagating in imaginary time. Fig. |^ 
shows examples of a 3D parent condensate wavefunction ^'(x, y, z, t) for two different times. 
The t = solution shows the wavefunction in the harmonic trap, and the t = ti = 600 
fis solution shows the wavefunction after 600 fis of free evolution without a trap potential. 
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Although the t = wavefunction in Fig. has a constant phase (taken to be 0), it is 
apparent from Fig. ^3 that the evolution leads to the development of phase modulation across 
the condensate, i. e., the wavefunction develops a spatially dependent phase, and therefore 
an imaginary part of the wavefunction. This is due to the evolution of the condensate under 
the influence of the mean field term, A^f/o|\l/(r, when the trapping potential is no longer 
present. An analytic form for the spatially dependent phase which evolves can be obtained 
in the Castin-Dum model |2^. As we show below, this phase modulation is important 
for 4WM. There is very little physical expansion of the condensate after 600 /is, since the 
condensate densities |\l/(r,t)p are nearly the same for the wavefunctions in Figs. ^ and ^b. 
However, Fig. ^ shows that the acceleration due to the mean field is already quite evident 
in the momentum distribution at t = 600 /is, which is much broader than that at t = 0. 
The two peaks near k = ±5r^p in the t = ti = 600 /xs distribution indicate the formation 
of accelerated condensate particles which will lead to condensate expansion at later times. 

Our treatment for applying Bragg pulses uses the model given by Eq. (^. This ap- 
proximation neglects detailed dynamics during the application of the Bragg pulses. Each 
initial wavepacket i at time t2 after the Bragg pulses is a copy of the parent condensate 
wavefunction at t = ti with population fraction = Ni/N. Unless stated otherwise, we 
will always use the ratio /i : /2 : /s = 7 : 3 : 7 of population fractions as typical of the 
NIST experiment 0. We let the three BEG wavepackets evolve for t > ^2 ~ using three 
different versions of the time-dependent GPE. Two of them are 2D versions, and one is the 
3D-SVEA version. The 2D-full version uses the GPE, Eq. (|l|), to evolve the initial state ^ 
in Eq. (^. The 2D-SVEA version uses the SVEA form in Eqs. (|T^)- ([T5|) for the evolution. 
A typical 2D calculation used a grid of discrete x, y points within a box 5rTF wide in the x 
and y directions centered on x = y = 0. In order to resolve the rapid phase variations due 
to the e*^'*^ '"-' factor, the 2D-full calculation required an x, y grid of up to 4096 x 4096 points. 
On the other hand, the 2D-SVEA only requires a 128 x 128 x, y grid to achieve comparable 
accuracy. The 3D-SVEA calculations added a 4rrF wide box in the z direction, and an 
X, y, z grid of 128 x 128 x 64 was sufficient. 

Fig. ^ compares the 4WM output fraction ji^if) = N4{t)/N for the three different types 
of calculation for the case of = 1.5 x 10^ atoms. The 2D-full and 2D-SVEA calculations 
give the same results within numerical accuracy and can not be distinguished on the graph. 
We take this to be a strong justification of the SVEA, and a strong indication that it will 
be equally trustworthy in the 3D calculations. In both 2D and 3D cases, the output grows 
quadratically at early time, as predicted by Eq. (JIB]). The arrows indicate the characteristic 
nonlinear time Inl and the collision time tcoi- In addition, the figure shows tcoi{x) = tcoi/V^- 
The latter is the time it takes wavepackets 1 and 2 to move so that they just touch at their 
Thomas-Fermi radii in the x direction. At that time wavepackets 1 and 2 no longer have 
significant overlap with each other, although they still have some overlap with wavepacket 
3. As the wavepackets begin to move apart, the output saturates near t — ^2 ~ tcoi{.x)/2 
and approaches its final value when t — ^2 ~ icoi- There is a significant difference between 
the 3D-SVEA and 2D-SVEA output fraction. The 4WM output is lower for the 3D case. 
This is because the nonlinear 4WM process depends on the spatial overlap of the moving 
wavepackets. The packets are not as well-overlapped geometrically in 3D as in the 2D model. 
Henceforth, all our calculations are 3D-SVEA ones, unless stated otherwise. 

Fig. H shows a sequence of contour images of the time evolution of the wavepackets from 
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the time the trap is turned off at t = to the time of separation of the four wavepackets. 
The contours show the 2;-integrated column density, J2i=i I ^i{x,y, z,t)\'^dz, from the 3D- 
SVEA calculation. (The constructive and destructive interference fringes in the wavepacket 
overlap region due to the e*'^ '" phase factors is not shown since it would require very high 
resolution to represent it with sufficient accuracy). Panel (a) shows the eigenstate density 
in the harmonic trap. Panel (b) shows the wavepacket at t = t2 just after the Bragg pulses 
have fired. Since there is negligible expansion in the density profile during the initial 600 /is 
of free evolution, the wavepacket is very similar to that in panel (a). However, we learned 
from Fig. ^ that a phase modulation has developed across the wavepacket. This does not 
show up in the density profile. Panel (c) for t — t2 = 190 fis indicates some initial motion 
by the moving wavepackets. In panel (d) the spread of the three wavepackets due to their 
different momenta is evident, and in panel (e) the separation of the 4WM wavepacket is 
clearly apparent. Panel (e) shows the four wavepackets after almost complete separation at 
t — t2 = 760 /is, which is larger than tcoi = 687/is. 

Fig. compares the output fraction N4{t)/N versus time for three different initial total 
atom numbers, = 0.2 x 10^, 1.5 x 10^ and 5.0 x 10^, and ti = 600 /xs. Again, at early 
times the quadratic dependence of the fraction as a function of time is clearly evident. After 
a quadratic rise at early time, the output saturates and even undergoes oscillations before 
finally settling down to a final value when t > tcoi- The oscillations of N4{t)/N in time 
develop and become more pronounced as the initial number of atoms increases. These are 
due to back-transfer from the i = 2 and 4 packets to the i = 1 and 3 packets due to the 
mutual coupling between the packets. A closer examination of the detailed time evolution 
shows that the transfer occurs on the trailing edge of the wavepackets where they are still 
substantially overlapped. When N is large enough, the wavepackets experience significant 
distortion in shape by the time they separate. The output fraction N4^(t)/N clearly increases 
with A^. 

Fig. ^ shows the output fraction N4{t)/N versus time for 1.5 x 10^ atoms for four different 
values of the free evolution time ti = ^s, 600 /xs, 1200 /is, and 1800 fis. The self-phase 
modulation resulting from the nonlinear self-energy interaction reduces the 4WM output 
as ti increases. This is analogous to the destruction of third harmonic generation due to 
self- and cross-phase modulation in nonlinear optics |2^, and occurs because the phase 
modulation destroys the phase matching that is necessary for 4WM to develop. For t > tcoi, 
the number of atoms in the different wavepackets no longer change, since the wavepackets 
are well separated (exchange of the number of bosonic atoms between wavepackets can no 
longer occur when the terms in the dynamical equations responsible for 4WM vanish). From 
these calculations it seems clear that 4WM should be much stronger if the trap is left on 
instead of being turned off. These calculations indicate that the 4WM output of the NIST 
experiment might be as much as a factor of two higher if there had not been 600 fis of 
free evolution before the Bragg pulses were applied. 

We expect the 4WM output will be larger if the wavepackets stay together for a longer 
interaction time tcoi- The interaction time can be changed by changing the velocity of the 
wavepackets. Fig. |^ plots Ni{t)/N versus time for 1.5 x 10^ atoms for the original case shown 
in Figs. 1^ and ^ and for two new cases where the interaction times are changed by factors 
of 0.7 and 2. This is achieved in the code by scaling the momentum wavevectors by factors 
of 1/0.7 and 1/2 respectively. Our calculations show that the 4WM output is reduced by 
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a factor of 0.6 in the first case and increased by a factor of 2 in the second. In principle, 
velocities of the wavepackets can be controlled by changing the frequencies and angle of the 
two Bragg pulses that create an outcoupled wavepacket 0. Thus, some degree of control 
over the 4WM output should be possible by varying the interaction time. 

Fig. |10| shows fsit) and f4{t) for the case of a weak i = 2 "probe" with initial population 
fraction 0.001 incident on two strong i = 1 and 3 "pump" wavepackets with population 
fractions 0.4995. This is analogous to the phase conjugation process envisioned in reference 
P]. Here bosonic stimulation, which removes 2 atoms from the "pump" packets 1 and 3 and 
puts them in packets 2 and 4, results in a strong amplification of packet 2, which grows in 
atom number 8-fold as the 4WM signal grows. 

Fig. |ll] shows 4WM output fraction N4/N after the half-collision is over (t > tcoi) as a 
function of A^, plotted in a log-log plot. The figure shows the results for both the 2D-SVEA 
and 3D-SVEA calculations. The dashed lines show the 4WM output for small N scales well 
with N^^^, as estimated from the simple model in Section [11 D\ The scaling with N^^^ for 
small N is clearly evident in both 2D and 3D results. The latter is uniformly lower than the 
former, due to the smaller overlap of the wavepackets in 3D because of geometrical reasons, 
but saturates a little more slowly with increasing than the former. At the higher N values 
typical of Na condensates, this scaling from the simple model seriously overestimates the 
output, which begins to saturate with increasing N. 

Fig. shows three curves giving the fraction of atoms in the 4WM output wavepacket 
as a function of the initial total number of atoms as calculated by (1) 2D-SVEA and (2) 
3D-SVEA simulations without including elastic scattering loss, and as calculated by (3) a 
3D-SVEA simulation including elastic scattering loss. In one set of calculations we used a 
ratio of atoms in the three initial wavepackets of A^i : A''2 : A'^a = 7 : 3 : 7. These calculations 
produce the three smooth curves in Figure |T2|. In another set of calculations, we used 
the measured final fractions from the NIST experiment |0| to determine the initial ratios 
Ni : N2 : N3, rather than taking the nominal values 7:3:7. The open circles in Figure 



12, which no longer fall on a smooth line, show the 3D-SVEA without elastic scattering for 



these cases with experimental scatter in initial conditions. The relatively small deviation of 
the points from the solid curve for the 3D-SVEA without elastic scattering show that the 
calculations with the 7:3:7 ratio is useful for generating a smooth curve to compare to 
experimental data. 

The effect of including loss from the BEC wavepackets due to elastic scattering collisions 
was modeled using Eqs. (pil|)-(P5|). The 4WM output reduction in Figure |T^ due to elastic 
scattering ranges from 6 per cent to 16 per cent in going from 10^ to 10^ atoms, and becomes 
more pronounced for large values of A^, with the loss due to elastic scattering reaching 36 per 
cent for 5 x 10® atoms. Elastic scattering of atoms from the different momentum wavepackets 
removes atoms from the four BEC wavepackets, and it thereby also lowers the nonlinear 
coupling term that gives rise to the 4WM. Although the mean-free-path for elastic collisions 
is on the order of 10 times ttf for 1.5 x 10® atoms, there are a sufficient number of collisions 
to make a noticable reduction in the nonlinear output. 

Finally, Fig. |13| compares our 3D-SVEA calculation, with corrections due to elastic scat- 
tering, to the observed output 4WM fraction in the NIST experiment . The overall agree- 
ment is good, given the approximations in the model and the scatter in the experimental 
data. The calculated curve tends to be slightly larger than the mean of the measured points, 
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and in particular, does not seem to saturate as fast at large as the experimental data. 
Since systematic error bars were not given for the data, it is difficult to know whether this 
slight disagreement is significant. There are clearly approximations in the theory, such as 
using the GPE method or ignoring the dynamics during the application of the Bragg pulses. 
There also are effects in the experiment that might have a bearing on the comparison. For 
example. Fig. 2b of reference reported a best case of 10.6 per cent 4WM output for 
= 1.7 X 10^ atoms, although a lower figure near 6 per cent reported in Fig. 3 of reference 
was more typical. The 10.6 per cent output would disagree with our calculations on the 
high side. This indicates that there is sufficient uncertainty in the quantitative aspects of 
the experiment to warrant a more systematic experimental exploration of the 4WM signal. 
Other possible sources of differences between theory and experiment include micromotion of 
the initial BEG in the time-orbiting-trap, laser misalignment, and a small finite temperature 
component of the BEG. 

IV. SUMMARY AND CONCLUSIONS AND OUTLOOK 

We have developed a full description of four- wave mixing (4WM) using a mean-field treat- 
ment of Bose- Einstein condensates. The slowly- varying-envelope approximation is a powerful 
tool that reduces the numerical grid requirements for calculating the time-dependent dynam- 
ics of fast-moving wavepackets with velocities greater than a photon recoil velocity. We find 
that elastic scattering loss between atoms in the fast wavepackets removes enough atoms 
from the wavepackets to affect the 4WM output. The quantum mechanical 3D calculations 
presented here show good agreement with experiment. 

In spite of the strong analogy between atom and optical 4WM, there are fundamental 
differences. In optical 4WM, the energy-momentum dispersion relation is different than 
in the massive boson case. Because we neither create nor destroy atoms, the only 4WM 
processes allowed for matter waves are particle number conserving. This is not the case for 
optical 4WM where, for example, in frequency tripling three photons are annihilated and 
one is created. Particle, energy and momentum conservation limit all matter 4WM processes 
to configurations that can be viewed as degenerate 4WM in an appropriate moving frame. 

We have considered 4WM using condensates of the same internal states. The internal 
states of the atoms can be changed by using Raman transitions. Thus, one can envision 
scattering atoms in one internal state from the matter-wave grating formed by atoms in a 
different internal hyperfine state. It is also possible to study the details of 4WM between 
mixed atomic species. We are in the process of carrying out such calculations. Quantum 
correlations created by the nonlinear process could lead to the study of non-classical matter- 
wave fields, analogous to squeezed and other non-classical states of light. It is of interest 
to investigate such cases. By varying the magnetic field to allow a Feshbach resonance to 
change the Uq coupling parameter, 4WM can be modified dynamically during the dynamics 
that occur as the wavepacket fiy apart, thus increasing or decreasing 4WM output. Such 
studies are also feasible. 

It is possible to modify the mean-field description of 4WM, and more generally, Bragg 
scattering of BEGs, by generalizing the GP equation to allow incorporation of momentum 
dependence of the nonlinear parameters, thereby putting the treatment of elastic and in- 
elastic scattering on a firm footing. This will be presented elsewhere | |26| ]. 
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FIG. 1. Momentum space view of the wavepackets participating in the four- wave mixing pro- 
cess, (a) Conservation of momentum in the laboratory frame, (b) A set of possible wavepackets 
in the laboratory frame with momenta that satisfy the phase-matching conditions in Section II C| , 
namely, IP2I = jPslcos^. 
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FIG. 2. (a) Lab frame view of the four-wave mixing process, showing the four wavepackets at 
early time while they are still interacting and at late time after they have separated, (b) Degenerate 
frame view of the same cases as in (a). 
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FIG. 3. (a) Cuts along the x, y and z axes of the parent condensate wavefunction 
^{x,y, z,t = 0) for N = 1.5 x 10^ atoms in a trap with harmonic frequencies of 84 Hz, 59.4 
Hz, and 42 Hz in the respective x, y, and z directions. The arrows show the TF radii rxpii) in the 
i = x,y,z directions. The curves labeled "x", "y", and "z" respectively represent Re[^'(x, 0, 0, 0)], 
Re[^(0,y,0,0)], and Rc[^(0, 0, z, 0)]; lm[*(x, y, z, 0)] is identically zero for each case, (b) Cuts 
along the x axis of Re[*(x, 0, 0, t = ti)] and Im[*(x, 0, 0, t = ti)] for ti = 600 /xs. 
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FIG. 4. Cut in the kx direction {ky = k^ = 0) of the squared momentum distribution |\I'(k, t)p 
for the wavefunctions in Fig. ^ for ti = and ti = 600 ^s. 
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FIG. 5. Comparison of N4(t)/N versus t — t2 for 2D and 3D calculations for 1.5 x 10® atoms. 
The trap is the same as in Fig. |3[ The Bragg pulses are applied 600 /is after the trapping potential 
is turned off and are over at time t2- 
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FIG. 6. Contour plots of integrated column density from the 3D-SVEA calculations vs x and y 

for iV = 1.5 xlO^ and the same trap as for Fig. |^. Panels (a) through (f ) show the time development 
of the wavepackets from the from the time the trap is turned off until the wavepackets physically 
separate. 



23 



0.20 



0.15 - 



IZ 0.10 



0.050 




0.0 

100 200 300 400 500 600 700 800 
t - \ (^s) 

FIG. 7. Comparison of Ni{t)/N versus t - ta for 0.2 x 10*^, 1.5 x 10*^ and 5.0 x 10^ atoms. The 
trap is the same as in Fig. |3|. The Bragg pulses are appUed 600 /xs after the trapping potential is 
turned off. 
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FIG. 8. Comparison of N4{t)/N versus t — t2 for 1.5 x 10^ atoms. The different curves show 
cases where the Bragg pulses are applied at ti = 0, 600, 1200 and 1800 fis after the trapping 
potential is turned off (ta ~ ^i)- The trap is the same as in Fig. |^. 
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FIG. 9. Comparison of N4{t)/N versus t — t2 for 1.5 x 10^ atoms. The trap is the same as in 
Fig. ^. The Bragg pulses are apphed 600 /xs after the trapping potential is turned off. The three 
different curves are for the cases where the separation times are scaled by factors of 0.7, 1, and 2 
by scaling the separation velocities by 1/0.7, 1, and 1/2. 
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FIG. 10. Growth of N4{t)/N and N4^{t)/N versus t — t2 for the case where a weak probe 
wavepacket 2 with initial population fraction 0.001 encounters strong "pump" wavepackets with 
initial fractions 0.4995. The trap is the same as in Fig. ^. The Bragg pulses are applied 600 ^s 
after the trapping potential is turned off. 
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FIG. 11. N/j^/N dependence dependence on the total number of atoms, A'^, calculated in 2D 
and 3D. The dashed lines show the N^/^ dependence predicted by the simple theory in subsection 
II D| . The trap is the same as in Fig. ^. The Bragg pulses are applied 600 fis after the trapping 
potential is turned off. 
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FIG. 12. Fraction of atoms in the 4WM output wavepacket, N4^/N, versus the total number of 
initial atoms, N, calculated in 2D, 3D and 3D with inclusion of elastic scattering loss as discussed 
in Sec. HE. The open circles represent calculations using experimental data to determine the 
ratios Ni : N2 : rather than taking the nominal values A'"! : N2 : = 7:3:7. The trap is the 
same as in Fig. 0. The Bragg pulses are applied 600 /is after the trapping potential is turned off. 
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FIG. 13. Fraction of atoms in the 4WM output wavepacket, N4^/N, versus the total number of 
initial atoms, N, calculated in 3D without and with inclusion of elastic scattering loss as discussed 
in Sec. IIE . The dots are experimental data Q. The trap is the same as in Fig. |. The Bragg 
pulses are applied 600 /.is after the trapping potential is turned off. 
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